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Abstract-  Most  current  computer-aided  detection  (CAD) 
algorithms  for  the  fully  automatic  detection  of  colonic  polyps 
from  3D  CT  data  suffer  from  high  false  positive  rates.  We 
developed  and  evaluated  a  post-processing  algorithm  to  decrease 
the  false  positive  rate  of  such  a  method.  Our  method  attempts  to 
model  the  way  a  radiologist  recognizes  a  polyp  while  scrolling  a 
cross-sectional  plane  through  3D  CT  data  by  quantifying  the 
change  in  location  of  the  edges  in  2D  plane.  It  uses  a  classifier 
for  identification  based  on  the  Mahalanobis  distance.  The  new 
method  increased  the  ROC  curve  area  from  0.89  to  0.98  (an 
increase  from  34.5%  to  85.0%  in  specificity  for  100% 
sensitivity)  in  a  population  of  8  patients. 

Keywords  -  Virtual  colonoscopy,  Computed  Tomographic 
Colonography,  Polyp  Detection,  CAD,  Optical  Flow  Fields 

I.  Introduction 

Computed  Tomographic  Colonography  (CTC)  is  a  minimally 
invasive  method  based  on  the  examination  of  the  colon  using 
CT  volume  data  [1].  This  has  been  achieved  by  examining 
the  stack  of  2D  images,  3D  virtual  colonoscopic  views,  or 
both.  Recently,  research  on  CTC  has  shifted  toward 
developing  fully  automatic  CAD  methods  for  polyp 
detection.  The  following  methods  have  been  proposed  for 
computer-aided  detection  (CAD)  of  colonic  polyps:  shape- 
based  detection,  Hough  Transform  (HT)-  based  detection  and 
sphere  fitting  [2-7].  Despite  their  high  sensitivity,  these 
methods  suffer  from  large  false  positive  rates  which  require 
radiologists  to  examine  a  large  number  of  images  manually. 
Reducing  the  number  of  false  positives  would  reduce  the 
radiologists’  reading  time  and  help  make  CTC  cost  effective. 

Our  goal  is  to  develop  and  validate  a  highly  sensitive  and 
specific  CAD  method,  intended  to  be  used  as  a  post-processor 
to  reduce  false  positives  from  a  high  sensitivity,  low 
specificity  polyp  detector.  We  attempted  to  model  the  way  a 
radiologist  recognizes  a  polyp  by  scrolling  through  a  stack  of 
2D  CT  images.  We  used  Optical  Flow  Fields  (OFFs)  to 
represent  changes  in  consecutive  images  and  a  linear 
classifier  to  sort  true  polyps  from  false  positives  based  on 
features  extracted  from  the  computed  OFFs. 

II.  Methodology 

A.  Pre-processing 

The  3D  CT  data  was  pre-processed  by  custom  software 
which  automatically  computes  a  path  along  the  central  colon 
axis  stretching  from  the  rectum  to  the  cecum  [8],  identifies 
the  colon  wall,  and  uses  a  HT  based  poylp  detector  (HT)  that 
detects  spherical  surface  patches  on  the  colon  wall  to  localize 
polyp-like  structures  [3].  For  computational  efficiency,  we 
extracted  sub  volumes  of  21x21x21  voxels  centered  on  the 


HT  detected  points  (HT_hits)  to  use  as  inputs  in  the  next 
stage. 

B.  Post-processing 


Post-processing  is  comprised  of  OFF  computation  to 
represent  the  changes  in  the  location  of  edges  in  the  CT 
images  (tissue/air  boundaries)  as  one  scrolls  through  the 
volume  CT  data,  and  the  morphological  characterization  of 
the  computed  OFFs.  Let  xy-plane  be  the  image  plane 
perpendicular  to  the  scrolling  direction,  z.  The  basic  OFF 
equation  for  a  constant  z  value,  is  [9]: 

VI-v  +  — =  0  (1) 

dz 

where  v(x,  y)  is  along  the  local  VI,  I(x,  y)  representing 
the  image.  As  the  scrolling  is  done  from  the  center  of  the 
sub  volume  (HT_hit)  outwards  in  both  the  -Z  and  +Z 
directions,  the  edges  of  spherical,  polyp-like  structures,  move 
inwards  on  the  image  plane.  This  consistency  in  their  motion 

is  required  as  Vz(x,y)’s  for  all  z  are  summed  and 


smoothed  to  get  a  single  OFF,  V D  (x,  y)  ,  associated  with  the 
current  sub  volume  and  scrolling  direction,  as  follows: 


vD(x,y)  =  Smooth 


Zvz(x’y) 


(2) 


V  Z 


The  smoothing  filter  used  is  a  3x3  rectangular  window  acting 
on  the  image  plane  pixelwise.  This  is  repeated  for 

De  {coronal,  sagittal, axial},  resulting  in  3  OFFs  that 
encode  information  in  3  orthogonal  scrolling  directions. 

A  parent  (Pn)  and  eight  child  nodes  (Cn’s)  are  determined 
on  each  OFF  for  its  characterization.  The  Pn  is  defined  to  be 


the  minimum  divergence  node  in  a  5x5  pixel  neighborhood 
of  the  HT_hit.  Cn’s  are  defined  to  be  the  points  5  pixels  away 
on  the  streamlines  incoming  to  immediate  neighbors  of  the 
Pn.  Fig.  1  shows  three  axial  CT  images  and  the  associated 


vay/al(x,  y)  •  The  local  Jacobian  matrix  of  VD(x,  y)  is 
used  to  compute  a  and  (3  parameters  defined  as  follows  [10]: 

r dv*  dvx  i 


c)x  dy 
dvy  dv 


dx  dy 


P  =  -traced)  ,  Q=  J  ,  (4) 


a- 


P  =  sign(P  -  40  J  P  -  4<2 


(5) 
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Axial  Image  #1 

[  + 


Axial  Image  #3 


Axial  Image  #2 

f  + 


Axial  Optical  Flow  Field 


Figure  1 .  Three  sequential  axial  images  around  a  HT_hit  (shown  by 
+)  and  the  associated  axial  optical  flow  field.  The  Pn  is  marked  with 
a  small  square,  the  Cn’s  are  marked  with  small  circles  on  the 
streamlines. 


The  parameters  a  and  (3  depend  on  the  local  properties  of 
the  OFF.  Their  ratio  uniquely  describes  the  local  topology. 
They  essentially  carry  the  information  in  the  eigenvalues  of 
the  characteristic  equation  of  J  . 

The  behavior  of  the  streamlines  is  quantified  by  using  two 
parameters  that  describe  the  spread  of  the  locations  of  Cn’s 
around  the  Pn.  The  first  parameter,  d  ,  is  defined  as 


,etj  =  Z(Cn1,CnJ.)e[0,;r]  (6). 


The  angle,  0- ,  is  calculated  with  respect  to  the  Pn’s  location. 

The  second  parameter,  Sc  ,  is  the  angle  of  the  largest  sector 
around  the  Pn  that  does  not  contain  a  Cn. 

Each  parameter  is  calculated  for  the  3  orthogonal  OFFs, 
making  up  a  12  dimensional  feature  vector. 


B.  Classification 


A  Mahalanobis  distance  based  classifier  is  used  for 
learning  and  classification  [11].  The  Mahalanobis  distance  of 

a  feature  vector  f  to  the  mean  vector,  mr ,  of  a  training  set, 
r ,  is  defined  as, 

rt,  mr  =  V(f -mr)TCr_I(f-mr)  (7) 

where  Cr  is  the  covariance  matrix  of  the  training  set.  This 
distance,  r,  is  a  standardized  measure  which  (i) 
automatically  accounts  for  scaling,  (ii)  takes  care  of 
correlations  between  features,  (iii)  can  provide  linear  and 
curved  decision  surfaces.  Referring  to  the  subset  of  polyps  in 
the  training  set  T  as  T1  and  the  subset  of  non-polyps  as 


Figure  2.  ROC  surface  for  a  single  experiment.  The  points  on  the  ROC 
surface  corresponding  to  \|/all  6|/ht)  are  marked  with  +  (*). 

ro  ,  the  feature  vector  fen  ,  where  n  is  the  test  set,  is 
classified  as 


f,mr 


-  rt 


f  ,mr 


+  b<  0 


otherwise 


fe  £11 
>f  e  £10 


(8) 


where  b  is  a  bias  term  and  £11  and  £10  are  polyp  and  non¬ 
polyp  subsets  of  £1  respectively. 


C.  Preliminary  Evaluation 


Data  were  acquired  from  8  patients  (7  male,  age  41-85, 
mean  age  63)  in  supine  position  after  colon  cleansing  and  air- 
insufflation.  Imaging  parameters  were  3mm  (2.5mm) 
collimation,  pitch  1. 5-2.0  (3.0),  1.5mm  (1.0- 1.5mm) 

reconstruction  interval,  200  MAs  (56  MAs),  and  120kVp,  on 
a  single-  (multi-)  detector  CT  system  (GE  Medical  Systems, 
Milwaukee,  WI).  Between  236  and  403  512x512  images 
were  reconstructed  for  each  case  with  an  average  voxel  size 
of  0.74mm  x  0.74mm  x  1.31mm.  These  patients  contained 
19  polyps  (5.0  to  23.0  cm.  diam.)  confirmed  by  fiberoptic 
colonoscopy. 

The  HT  detector  produces  a  score  related  to  the  size  and 
sphericity  of  each  HT_hit.  We  excluded  all  HT_hits  with  a 
score  below  an  arbitrary  threshold  h  =  10000,  resulting  in 
219  HT_hits  (19  true  polyps,  200  false  positives).  The 
sub  volumes  (21x21x21  voxels)  around  each  HT_hit  were 
used  as  the  input  data  to  the  OFF-based  classifier. 

The  preprocessed  data  set  was  divided  into  10  subsets 
randomly  such  that  the  true  positive  and  the  false  positive 
HT_hits  were  distributed  uniformly.  We  conducted  a  10-fold 
cross  validation  study  in  which  one  of  these  subsets  was  used 
as  the  test  set  and  the  rest  as  the  training  set  for  each 
experiment. 

The  sensitivity  and  specificity  are  defined  as  the 
percentages  of  the  correctly  detected  polyps  and  non-polyps, 
respectively,  within  the  test  set.  A  multidimensional  ROC 
surface  is  constructed  for  each  experiment  as  follows:  The 
HT_hits  corresponding  to  a  threshold,  h ,  are  fed  into  the 
OFF-based  classifier  (OFFC)  and  an  ROC  curve  is  computed 
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F# 

Param. 

Mean 

Median 

MinMax 

Sorted 

LeftMost 

Used 

A 

B 

C 

D 

E 

HT 

0.89 

fl 

[otpd] 

0.98 

0.97 

0.96 

0.96 

0.97 

f2 

[a  P] 

0.98 

0.97 

0.96 

0.97 

0.97 

f3 

[ad] 

0.97 

0.97 

0.97 

0.96 

0.97 

f4 

[pd] 

0.96 

0.95 

0.95 

0.95 

0.94 

f5 

[a  p  Sc] 

0.97 

0.97 

0.97 

0.96 

0.97 

f6 

[a  Sc] 

0.98 

0.97 

0.97 

0.97 

0.97 

f7 

[p  Sc] 

0.95 

0.96 

0.95 

0.96 

0.94 

Table  1.  ROC  areas  for  different  feature  sets:  i)  ‘Mean’  and  ‘Median’ 
mean  that  the  average  or  median  value  of  each  parameter  among  3  directions 
is  used,  ii)  ‘MinMax’  means  that  the  minimum(a),  minimum(lpl), 
maximum(d)  and  minimum(Sc)  are  used,  iii)  ‘Sorted’  means  that  all 
parameters  for  all  3  directions  are  used  after  sorting,  i.e.  The  feature  vectors 
are  9D  or  6D  depending  on  the  number  of  parameters  used,  iv)  ‘LeftMost’ 
means  that  the  scrolling  direction  with  a  leftmost  (a,p)  point  (minimum 
cos(a  /  ll[a  p]ll))  on  ap-plane  is  used. 

by  varying  b .  This  is  repeated  for  all  h , 
h  E  [10000,  max_  HT  _  score \ .  Varying  h  relates  SnS 
to  SpS_HT,  where  SnS  and  SpS_HT  represent  the  sensitivity 
and  specificity,  respectively,  of  the  HT  detector  at  each  h  . 
Now  consider  a  given  SpS_HT  obtained  for  some  threshold 
h  .  Varying  b  results  in  a  new  ROC  curve  relating  SnS  to 
SpS_ALL.  It  represents  the  sensitivity  and  specificity  of  the 
combination  HT  followed  by  OFFC.  SnS  on  this  new  ROC 
curve  represents  the  sensitivity  of  that  combination  and 
SpS_ALL  represents  its  specificity.  Plotting  these  these 
curves  for  every  SpS_HT  constructs  an  ROC  surface  in  3D 
that  depicts  SnS  as  a  function  of  SpS_HT  and  SpS_ALL. 
Note  that  the  surface  is  not  defined  for  SpS_ALL  less  than 
SpS_HT  because  the  OFFC  cannot  increase  the  number  of 
false  positives  provided  by  the  preceding  HT  stage.  The  SnS 
versus  the  SPS_HT  as  one  reads  while  walking  along  the 
(SpS_AFF=SpS_HT)  line  (\|/HT)  gives  the  ROC  curve  of  HT 
alone  (ROC_HT).  If  the  set  of  points  VALL^re  defined  as, 

^ ALL  ~  {  SPS1  ’  SPS2  | 

SpS2  =  maxSpS  ALL(snS(SpS1 ,  SpS , )  =  SnS(SpSj ,  SpS2)), 
VSpS, }  (9) 

then  these  points  correspond  to  the  maximum  specificity  one 
can  get  by  post-processing  the  HT_hits(  h ),  without 
sacrificing  SnS.  In  other  words,  SnS  versus  the  SpS_AFF  as 
one  reads  while  walking  along  \|/all  gives  the  ROC  curve  of 
the  overall  system  (ROC_AFF).  Further  more,  the  area 
between  \|/all  and  \\fm  depicts  the  improvement  in  specificity 
due  to  OFFC.  Fig.  2  shows  a  ROC  surface  and  the  associated 
Vht  and  \|/all- 

ROC_HTs  and  ROC_AFFs,  thus  computed  for  each  of  10 
experiments,  are  averaged  in  both  sensitivity  and  specificity 
to  get  mROC_HT  and  mROC_AFF.  Note  that  every  point  on 
an  ROC_AFF  corresponds  to  a  ( h  ,  b  )  pair.  However,  while 
the  corresponding  points  on  different  ROC_AFFs  have  the 
same  h  ,  this  is  not  necessarily  true  for  b  ,  because  for  each 
point  of  an  ROC_AFF,  b  is  selected  such  that  the  specificity 


P-Values  0  34  0  05  q  qq 


Figure  3.  The  final  ROC  curves  of  the  HT  alone  (mROC_HT)  and 
the  full  system  (mROC_ALL),  using  flA,  averaged  over  10  cross- 
validation  experiments.  The  error  bars  indicate  the  variation  among 
10  experiments.  The  p- values  indicate  the  significance  of  difference 
between  the  two  curves  (p-value<0.05  is  taken  to  be  statistically 
significant) 

is  maximum  without  sacrificing  the  sensitivity.  So,  despite 
the  fact  that  the  term  “ROC  curve”  is  not  totally  valid  for 
mROC_AFF,  it  has  similar  properties  in  terms  of  allowing 
the  comparison  of  sensitivity  and  specificity.  While  the  true 
detection  performance  is  represented  by  an  average  ROC 
surface,  this  does  not  allow  us  to  assess  the  improvement  in 
specificity,  as  further  discussed  in  Section  IV.  Our 
assessment  is  based  on  the  comparison  of  these  curves  for 
different  feature  vectors  listed  in  Table  1  and  the  curve  for 
HT  alone.  We  compared  the  ROC  areas  and  the  results  of  the 
t-tests  performed  for  each  point  on  the  mean  ROC  curves 
separately,  with  a  null  hypothesis  that  the  mean  of  the 
specificity  differences  between  two  mean  ROC  curves 
(mROC_HT  and  mROC_AFF)  for  the  corresponding 
sensitivity  is  zero  over  10  cross-validation  experiments.  Thus 
the  p-value  for  a  certain  sensitivity  value  assesses  the 
significance  of  the  difference  between  specificity  values  of 
two  methods  at  that  sensitivity. 

III.  RESUFTS 

Table  1  summarizes  the  ROC  areas  obtained  for  different 
feature  vectors  and  the  HT  alone.  The  ROC  area  for  the  HT 
alone  was  0.89.  There  is  not  much  difference  between  the 
ROC  areas  for  different  choices  of  feature  vectors  (Mean 
ROC  area  =  0.96±0.01).  However,  there  is  an  average 
difference  of  0.07  between  the  areas  of  these  curves  and  that 
of  the  HT  alone.  The  point-wise  t-tests  between  all  feature 
sets  and  the  HT  alone  resulted  in  mean  p-values  (averaged 
over  the  whole  spectrum  of  sensitivity  values)  less  than  0.009 
(maximum  mean  p-value=0.009  was  between  HT  alone  and 
flC).  This  means  that  on  the  average  the  null  hypothesis  (i.e. 
the  mean  difference  between  the  specificities  of  the  ROC 
curves  is  zero)  can  be  rejected  for  each  sensitivity  value. 
Actually,  a  closer  look  at  these  p-values  for  individual 
comparisons  show  that  all  p-values  are  less  than  0.05  for 
sensitivities  greater  than  or  equal  to  0.25.  Fig.  3  shows  the 
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Figure  4.  The  average  ROC  surface,  computed  over  10  cross-validation 
experiments.  The  volume  under  the  ROC  surface  is  0.82 

result  of  one  such  comparison  performed  between  the  HT 
alone  and  flA. 

IV.  DISCUSSION 

These  results  suggest  that  the  ROC  area  increases  due  to 
the  OFF-based  post-processing  algorithm  and  that  the 
difference  between  the  ROC  curves  before  and  after  post¬ 
processing  is  statistically  significant.  There  is  a  consistent 
decrease  in  the  ROC  area  when  the  parameter  a  is  excluded 
(except  for  f7C),  suggesting  that  a  carries  the  most 
information.  The  parameters  d  and  Sc  do  not  seem  to  carry 
much  information. 

Based  on  these  results,  flA  seems  to  be  the  best  choice 
because  it  has  a  large  ROC  area,  and  it  uses  three  parameters 
designed  to  represent  different  qualities  of  the  OFFs.  Fig.  3 
shows  the  mROC_HT  and  mROC_AFF  using  flA,  together 
with  the  p-values  computed  for  every  point  on  ROC  curves, 
over  the  set  of  10  experiments.  For  flA,  the  OFF  classifier 
increased  the  ROC  area  from  0.89  to  0.98.  At  a  sensitivity 
level  of  100%  (95%)  the  specificity  increased  from  34.5% 
(75.5%)  to  85%  (93%),  which  corresponds  to  a  decrease  in 
the  average  number  of  false  positives  per  test  set  from  13.1 
(4.9)  to  3.0  (1.4). 

As  mentioned  at  the  end  of  Section  II,  mROC_AFF  is  not 
a  conventional  ROC  curve.  It  represents  the  average 
improvement  in  specificity  without  decreasing  sensitivity.  It 
should  be  interpreted  as  the  average  of  the  maximum 
specificity  achievable  by  the  full  system  for  a  given 
sensitivity.  The  performance  of  the  overall  system  is  actually 
given  by  the  ROC  surfaces  as  in  Fig.  2.  The  average  of  these 
surfaces  for  flA,  computed  in  all  three  dimensions  point- 
wise,  is  given  in  Fig.  4.  The  ROC  volume,  as  a  percentage  of 
the  maximum  possible  volume  (which  is  0.5  as  the  ROC 
surface  is  defined  in  half  of  the  whole  volume)  is  0.82. 

V.  CONCFUSION 

Our  method  for  computer-aided  detection  of  colonic 
polyps  was  inspired  by  the  processes  radiologists  go  through 


while  investigating  CT  colonography  data.  Namely,  they 
scroll  through  the  3D  CT  data  in  a  certain  direction(s), 
usually  only  perpendicular  to  the  axial  plane  but  also 
perpendicular  to  the  coronal  and  sagittal  planes  when 
necessary,  following  the  edge  information  in  2D  images.  In 
our  algorithm,  this  information  is  represented  as  an  OFF 
vector  field  and  the  resultant  OFF  is  characterized  by 
topology  dependent  parameters.  We  showed  that  our  post¬ 
processing  algorithm  improves  the  specificity  of  a  high 
sensitivity,  low  specificity  CAD  without  sacrificing 
sensitivity.  The  design  of  OFFC  does  not  assume  anything 
about  the  nature  of  the  pre-processing,  and  thus  we  expect  it 
to  be  applicable  to  any  high  sensitivity,  low  specificity  CAD 
output.  Our  results  show  that  there  is  relevant  information  in 
the  motion  of  edges  in  terms  of  polyp  detection  in  CTC. 
Further  research  is  required  to  validate  the  presented  method 
and/or  develop  new  ways  of  representing  this  information. 
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